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Abstract 

We derive the mean- field equations arising as the limit of a network of interacting 
spiking neurons, as the number of neurons goes to infinity. The neurons belong 
to a fixed number of populations and are represented either by the Hodgkin- 
Huxley model or by one of its simplified version, the Fitzhugh-Nagumo model. 
The synapses between neurons are either electrical or chemical. The network 
is assumed to be fully connected. The maximum conductances vary randomly. 
Under the condition that all neurons initial conditions are drawn independently 
from the same law that depends only on the population they belong to, we 
prove that a propagation of chaos phenomenon takes places, namely that in 
the mean- field limit, any finite number of neurons become independent and, 
within each population, have the same probability distribution. This proba- 
bility distribution is solution of a set of implicit equations, either nonlinear 
stochastic differential equations resembling the McKean-Vlasov equations, or 
non-local partial differential equations resembling the McKean-Vlasov-Fokker- 
Planck equations. We prove the well-posedness of these equations, i.e. the 
existence and uniqueness of a solution. We also show the results of some pre- 
liminary numerical experiments that indicate that the mean-field equations are 
a good representation of the mean activity of a finite size network, even for 
modest sizes. These experiment also indicate that the McKean-Vlasov-Fokker- 
Planck equations may be a good way to understand the mean-field dynamics 
through, e.g., a bifurcation analysis. 

* email: f irstname .nameOinria.f r 

'''College de France, Center of Interdisciplinary Research in Biology,, CNRS UMR 7241m 
INSERM U1050, Universite Pierre et Marie Curie ED 158. MEMOLIFE Laboratory of excel- 
lence and Paris Science Lettre. 11, place Marcelin Berthelot, 75005 Paris, France. 



AMS 2010 subject classification: 60F99, 60B10, 92B20, 82C32, 82C80, 
35Q80 



Keywords and phrases: Mean-Field Limits, Propagation of chaos, Stochas- 
tic differential equations, McKean-Vlasov equations, Fokker-Planck equations. 
Neural Networks, Neural assemblies. Ho dgkin- Huxley neurons, Fitzhugh-Nagumo 
neurons. 

1 Introduction 

Cortical activity displays highly complex behaviors which are often character- 
ized by the presence of noise. Reliable responses to specific stimuli often arise at 
the level of population assembhes (cortical areas or cortical columns) featuring a 
very large number of neuronal cells, each of these presenting a highly nonlinear 
behavior, that are interconnected in a very intricate fashion. Understanding the 
global behavior of large-scale neural assemblies has been a great endeavor in the 
past decades. One of the main interests of large scale modeling is character- 
izing brain function, which most imaging techniques are recording. Moreover, 
anatomical data recorded in the cortex reveal the existence of structures, such 
as the cortical columns, with a diameter of about 50/im to 1mm, containing 
of the order of one hundred to one hundred thousand neurons belonging to a 
few different types. These columns have specific functions. For example, in the 
human visual area VI, they respond to preferential orientations of bar-shaped 
visual stimuli. In this case, information processing does not occur at the scale of 
individual neurons but rather corresponds to an activity integrating the individ- 
ual dynamics of many interacting neurons and resulting in a mesoscopic signal 
arising through averaging effects, and that effectively depends on a few effective 
control parameters. This vision, inherited from statistical physics requires that 
the space scale be large enough to include sufficiently many neurons and small 
enough so that the region considered is homogeneous. This is in effect the case 
of cortical columns. 

In the field of mathematics, studying the limits of systems of particles sys- 
tems in interaction has been a longstanding problem and presents many tech- 
nical difficulties. One of the questions addressed in mathematics was to char- 
acterize the limit of the probability distribution of an infinite set of interacting 
diffusion processes, and the fluctuations around the limit for a finite number of 
processes. The first breakthroughs to find answers to this question are due to 
Henry McKean (see e.g. [53[ |52]). It was then investigated in various con- 
texts by a large number of authors such as Braun-Hepp [10], Dawson |23J, 
Dobrushin [25 , and most of the theory was achieved by Tanaka and collab- 
orators [70 , 69J7ll[72], and of course Alain-Sol Sznitman gll |Ml EH [66] . When 
considering that all particles (in our case neurons) have the same, independent 
initial condition, they mathematically proved using stochastic theory (Wasser- 
stein distance, large deviation techniques) that in the limit where the number 
of particles tend to infinity, any finite number of particles behave independently 
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of the other ones, and they ah present the same probabihty distribution, which 
satisfies a nonlinear Markov equation. Finite-size fluctuations around the limit 
are derived in a general case in [65] . Most of these models use standard hypoth- 
esis of global Lipschitz continuity and linear growth condition of the drift and 
diffusion coefficients of the diffusions, as well as Lipschitz continuity of the inter- 
action function. Extensions to discontinuous cadlag processes including singular 
interactions (through a local time process) were developed in [66] . Problems in- 
volving singular interaction variables (e.g. nonsmooth functions) are also widely 
studied in the field, but are not relevant in our case. 

In the present article, we apply this mathematical approach to the problem 
of interacting neurons arising in neuroscience. To this end we extend the theory 
to encompass a wider class of models. This implies the use of locally (instead of 
globally) Lipschitz coefficients and of a Lyapunov-like growth condition replac- 
ing the customary linear growth assumption for some of the functions appearing 
in the equations. The contributions of this article are threefold: 

1. We derive in a rigorous manner the mean- field equations resulting from 
the interaction of infinitely many neurons in the case of widely accepted 
models of spiking neurons and synapses. 

2. We prove a propagation of chaos property which shows that in the mean- 
field limit neurons become independent, in agreement with some recent 
experimental work [26 and with the idea that the brain processes infor- 
mation in a somewhat optimal way. 

3. We show numerically that the mean- field limit is a good approximation 
of the mean activity of the network even for fairly small sizes of neuronal 
populations. 

4. We suggest numerically that the changes in the dynamics of the mean- 
field limit when varying parameters can be understood by studying the 
mean-field Fokker-Planck equation. 

We start by reviewing such models in section [2] to motivate the present study. 
It is in section [3] that we provide the limit equations describing the behaviours 
of an infinite number of interacting neurons, and state and prove the existence 
and uniqueness of solutions in the case of conductance-based models. The de- 
tailed proof of the second main theorem, that of the convergence of the network 
equations to the mean- field limit, is given in appendix [A] In section [4| we begin 
to address the difficult problem of the numerical simulation of the mean-field 
equations and show some results indicating that they may be an efficient way 
of representing the mean activity of a finite-size network as well as to study the 
changes in the dynamics when varying biological parameters. The final section [5] 
focuses on the conclusions of our mathematical and numerical results and raises 
some important questions for future work. 
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2 Spiking Conductance-based models 



This section sets the stage for our results. We review in section 2.1 the Hodgkin- 
Huxley model equations in the case where both the membrane potential and the 
ion channels equations include noise. We then proceed in section 2.2 with the 
Fitzhugh-Nagumo equations in the case where the membrane potential equa- 
tion includes noise. We next discuss in section 2.4 the connectivity models of 
networks of such neurons, starting with the synapses, electrical and chemical, 
and finishing with several stochastic models of the synaptic weights. In section 



2.5 we write the network equations in the various cases considered in the pre- 



vious section and express them in a general abstract mathematical form that 
is the one used for stating and proving the results about the mean-field limits 
in section |3l Before we jump into this we conclude in section [2^ with a brief 
overview of the mean-field methods popular in computational neuroscience. 

From the mathematical point of view, each neuron is a complex system, 
whose dynamics is often described by a set of stochastic nonlinear differen- 
tial equations. Such models aim at reproducing the biophysics of ion channels 
governing the membrane potential and therefore the spike emission. This is 
the case of the classical model of Hodgkin and Huxley [41 and of its reduc- 
tions [35l |36l [42] . Simpler models use discontinuous processes mimicking the 
spike emission by modeling the membrane voltage and considering that spikes 
are emitted when it reaches a given threshold. These are called integrate-and-fire 
models ^47^ i78j and will not be addressed here. The models of large networks 
we deal with here therefore consist of systems of coupled nonlinear diffusion 
processes. 



2.1 Hodgkin-Huxley model 

One of the most important models in computational neuroscience is the Hodgkin- 
Huxley model. Using pioneering experimental techniques of that time, Hodgkin 
and Huxley [4r determined that the activity of the giant squid axon is controlled 
by three major currents: voltage-gated persistent current with four activa- 
tion gates, voltage-gated transient Na^ current with three activation gates and 
one inactivation gate, and Ohmic leak current, /l, which is carried mostly by 
chloride ions Cl~ . In this paper we only use the space-clamped Hodgkin-Huxley 
model which we slightly generalize to a stochastic setting in order to better take 
into account the variability of the parameters. 

The basic electrical relation between the membrane potential and the cur- 
rents is simply: 

= ^'"'(i) -Ik- iNa - II, 
jext^^j is an external current. The detailed expressions for i^, I^a and II can 
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be found in several textbooks, e.g., [42| l3Q]: 

lK=9Kn\V -Ek) 

iNa = QNam^hiy - Eno) 

Il=9l{V-El) 

where qk (respectively Qno) is the maximum conductance of the potassium 
(respectively the sodium) channel, ql is the conductance of the Ohmic chan- 
nel, n (respectively m) is the activation variable for (respectively for No). 
There are four (respectively 3) activation gates for the (respectively the 
No) current which accounts for the power 4 (respectively 3) in the expression 
of Ik (respectively In a)- h is the inactivation variable for Na. These activa- 
tion/deactivation variables, denoted by x G {n, m, h} in what follows represent 
a proportion of open channels (they vary between and 1). This proportion 
can be computed through a Markov chain modeling assuming the channels to 
open with rate Px{V) (the dependence in V accounts for the voltage-gating of 
the channel) and to close with rate Cx{V). These processes can be shown to con- 
verge, under standard assumptions, towards the following ordinary differential 
equations: 

X = Px{V){l - x) - Cx{V) X, X e {n, m, h} 

The functions px{V) and Cx{V) are smooth and bounded functions whose exact 
values can be found in several textbooks such as the ones cited above. A more 
precise model taking into account the finite number of channels results in the 
stochastic differential equation (see e.g. [81 and references therein): 

dxt = {Px{V){l -x)- UV) x) dt + ^/px{V){l-x)^UV)x x{x) dWr 

where I^f , x G {n, m, h} are independent standard Brownian motions, x(^) 
is a function that vanishes outside (0,1). This guarantees that the solution 
remains a proportion i.e. lies between and 1 for all times. We define 

a,{V,x) = y'p,{V){l-x) + UV)x x{x). (1) 



In order to complete our stochastic Hodgkin-Huxley neuron model, we as- 
sume that the external current P^^{t) is the sum of a deterministic part, noted 
/(t), and a stochastic part, a white noise with variance ciext built from a stan- 
dard Brownian motion Wt independent of Wf , x G {n, m, h}. Considering the 
current produced by the income of ion through these channels, we end up with 
the following stochastic differential equation: 

^CdVt = (/(t) - SKU^iV - Ek) - SNam^hiV - Eno) - gL{V - El)) dt 

< ^CFext dWt 

dxt = (^Px{V){l - x) - Cx{V)x) dt + ax{V,x) dW^ x G {n, m, h} 

(2) 
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The advantages of this model are numerous, and one of the most prominent 
aspects in its favor is its correspondence with the most widely accepted for- 
malism to describe the dynamics of the nerve cell membrane. A very extensive 
literature can also be found about the mathematical properties of this system, 
and it is now quite well understood. 

The stochastic version of the space-clamped Hodgkin-Huxley model we study 
is, from the mathematical point of view, a four-dimensional nonlinear stochastic 
differential equation. It is easy to check that the drift and diffusion functions are 
all Lipchitz-continuous and hence satisfy the linear growth condition, using the 
fact that all variables n, m and h are in the interval [0, 1] (we recall that these 
variables are proportions of open channels). Standard theorems of stochastic 
calculus ensure under these two conditions existence and uniqueness of strong 
solutions for any sufficiently regular (square integrable) initial condition (see 
e.g. [45, 50 ). 

To illustrate the model we show in figure [l] the time evolution of the three 
ion channel variables n, m and h as well as that of the membrane potential V 
for a constant input / = 20.0. The system of ODEs has been solved using a 
Runge-Kutta 4 scheme with integration time step At = 0.01. 

In figure |2] we show the same time evolution when noise is added to the 
channel variables and the membrane potential. 

For the membrane potential, we have used a^xt = 3.0 (see equation ([2|), while 
for the noise in the ion channels we have used the following x function (see 
equation ([T]) 
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if < X < 1 
ifx<OVx>l 



(3) 



with r = 0.1 and A = 0.5 for all the ion channels. The system of SDEs has 
been integrated using the Euler-Maruyama scheme with At = 0.01. 



Evolution of a single Hodgkin-Huxley neuron without noise 
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Figure 1: Solution of the noiseless Hodgkin-Huxley model. Left: Time evolution 
of the three ion channel variables n, m and h. Right: Corresponding time 
evolution of the membrane potential. Parameters given in the text. 
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Figure 2: Noisy Ho dgkin- Huxley model. Left: Time evolution of the three 
ion channel variables n, m and h. Right: Corresponding time evolution of the 
membrane potential. Parameters given in the text. 



Because the Hodgkin-Huxley model is rather complicated and high-dimensional, 
many reductions have been proposed, in particular to two dimensions instead of 
four. These reduced models include the famous Fitzhugh-Nagumo and Morris- 
Lecar models. These two models are two-dimensional approximations of the 
original Hodgkin-Huxley model based on quantitative observations of the time 
scale of the dynamics of each variable and identification of variables. Most 
reduced models still comply with the Lipschitz and linear growth conditions en- 
suring existence and uniqueness of a solution, except for the Fitzhugh-Nagumo 
model which we now introduce. 



2.2 The Fitzhugh-Nagumo model 

In order to reduce the dimension of the Hodgkin-Huxley model, Richard Fitzhugh 
j34| \35\ [36] introduced a simplified two-dimensional model. The motivation was 
to isolate conceptually essential mathematical features yielding excitation and 
transmission properties from the analysis of the biophysics of sodium and potas- 
sium fiows. Nagumo and collaborators [55^ followed up with an electrical system 
reproducing the dynamics of this model and studied its properties. The model 
consists of two equations, one governing a voltage-like variable V having a cubic 
nonlinearity and a slower recovery variable w. It can be written as: 

V = f{V) -w + 

w = c{V + a — bw) 

where f{V) is a cubic polynomial in V which we choose, w.l.o.g., to be f{V) = 
V — The parameter models the input current the neuron receives, 

the parameters a, 6 > and c > describe the kinetics of the recovery variable 
w. As in the case of the Hodgkin-Huxley model, the current is assumed 
to be the sum of a deterministic part, noted /, and a stochastic white noise 



7 



accounting for the randomness of the environment. The stochastic Fitzhugh- 
Nagumo equation is deduced from Q and reads: 



\dVt = {Vt-^-Wt^I)dt^a,^tdWt 
]dwt = c{Vt -\- a — bwt) dt 

Because the function f{V) is not globally Lipschitz continuous (only locahy), 
we prove the fohowing proposition stating that the stochastic differential equa- 
tion ([5| is well-posed. 

Proposition 1. LetT > be a fixed time. If\I{t)\ < Im on [0,T]^ equation ^ 
together with an initial condition (Vo^wq) with a given distribution has a unique 
strong solution which belongs to L^([0, T]; R^). 

Proof. The function h = (/ii,/i2) : {V,w) {f{V) - w ^ Id, c{V ^ a - bw)) 
is locally Lipschitz in R^. According to [50, Chapter 2, Theorem 3.5] it is 
sufficient to show that, if we note X the vector (V, i^) of R^, there exists a 
positive constant K such that for all (X, t) in R^ x [0,T] 

X^h{X,t) + \al^,<K(l + \X\^) 
Let L = (c — l)/2, we have: 

X^h{X,t) = V^-X^+{c-l)V w-bcw'^+V I+acw < V^+2LVw+2V^+2^w 
<{V^Lwf^(v^^-]\{w^ac)^<2(2V^^{L^^l)w^ ' + a V 



2 J ' ' - \ ' ' 4 

r2 



< 2 max |2, + 1, ^!^!^ | (1 + F 



2 I .,,2 



The conclusion follows since / is bounded. □ 

We show in figure [S] the time evolution of the adaptation variable and the mem- 
brane potential in the case where the input / is constant and equal to 0.7. The 
lefthand side of the figure shows the case with no noise while the righthand side 
shows the case where noise of intensity ciext = 0.25 (see equation ([5|) has been 
added. 

The deterministic model has been solved with a Runge-Kutta 4 method, while 
the stochastic model with the Euler-Maruyama scheme. In both cases, we have 
used an integration time step At = 0.01. 



2.3 Partial conclusion 

We have reviewed two main models of space-clamped single neurons, the Hodgkin- 
Huxley and the Fitzhugh-Nagumo models. These models are stochastic, in- 
cluding various sources of noise, external, and internal. The noise sources are 
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Evolution of a single FitzHugln-rJagumo neuron witli noise 
P Evolution of a single FitzHugh-Nagumo neuron without noise 3\ ! 1 < < 




Time Time 



Figure 3: Time evolution of the membrane potential and the adaptation variable 
in the Fitzhugh-Nagumo model. Left: without noise. Right: with noise. See 
text. 

supposed to be independent Brownian processes. We have shown that the re- 
sulting stochastic differential equations ([2| and ([5| were well-posed. As pointed 
out above, this analysis extends to a large number of reduced versions of the 
Hodgkin-Huxley such as those that can be found in the book [42 . 

2.4 Models of synapses and maximum conductances 

We now study the situation in which several of these neurons are connected to 
one another forming a network, which we will assume to be fully connected. 
Let N be the total number of neurons. These neurons belong to P populations, 
e.g. pyramidal cells, interneurons. ... If the index of a neuron is z, 1 < z < A^, 
we note p{i) = a, 1 < a < P the population it belongs to. We note A'p(i) 
the number of neurons in population p{i). Since we want to be as close to 
biology as possible while keeping the possibility of a mathematical analysis of 
the resulting model, we consider two types of simplified, but realistic, synapses, 
chemical and electrical or gap junctions. The following material concerning 
synapses is standard and can be found in textbooks [30]. The new, and we 
think important, twist is add noise to our models. To unify notations, in what 
follows i is the index of a postsynaptic neuron belonging to population a = p{i)^ 
and j is the index of a presynaptic neuron to neuron i belonging to population 
7 =P(J)• 
2.4.1 Chemical synapses 

The principle of functioning of chemical synapses is based on the release of 
neurotransmitter in the presynaptic neuron synaptic button, which binds to 
specific receptors on the postsynaptic cell. This process is, similarly to the 
currents described in the Hodgkin and Huxley model, governed by the value 
of the cell membrane potential. We use the model described in ^241 [30], which 
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features a quite realistic biophysical representation of the processes at work in 
the spike transmission, and is consistent with the previous formalism used to 
describe the conductances of other ion channels. The model emulates the fact 
that following the arrival of an action potential at the presynaptic terminal, 
neurotransmitter is released in the synaptic cleft and binds to postsynaptic 
receptor with a first order kinetic scheme. Let j be a presynaptic neuron to the 
postynaptic neuron i. The synaptic current induced by the synapse from j to i 
can be modeled as the product of a conductance gij with a voltage difference: 

I:^=9^mv'-vL) (6) 

The synaptic reversal potential F/g^ is approximately constant within each pop- 
ulation: V^\^ := y^g^. The conductance gij is the product of the maximum 
conductance Jij {t) with a function (t) that denotes the fraction of open chan- 
nels and depends only upon the presynaptic neuron j: 

9ijit) = JijWvHt) (7) 

The function (t) is often modelled [30] as satisfying the following ordinary 
differential equation 

yHt)=aiS,{V^){l-yHt))-aiyHt), 

The positive constants and a;^ characterize the rise and decay rates, respec- 
tively, of the synaptic conductance. Their values depend only on the population 
of the presynaptic neuron j, i.e. := and a;^ := a^, but may vary signif- 
icantly from one population to the next. For example, GABAb synapses are 
slow to activate and slow to turn off while the reverse is true for GABAa and 
AMPA synapses [30 . Sj{V^) denotes the concentration of transmitter released 
into the synaptic cleft by a presynaptic spike. We assume that the function Sj 
is sigmoidal and that its exact form depends only upon the population of the 
neuron j. Its expression is given by (see e.g. [30 ): 

Destexhe et al. [24 give some typical values of the parameters, T^ax = ImM^ 
Vt = 2mV, and 1/A = 5mV. 

Because of the dynamics of ion channels and of their finite number, similarly 
to the channel noise models derived through the Langevin approximation in the 
Hodgkin-Huxley model ([2|, we assume that the proportion of active channels is 
actually governed by a stochastic differential equation with diffusion coefficient 
ary{V,y) depending only on the population 7 of j of the form ([T]): 

dyi = {aJS^{V^){l - y^{t)) - {t)) dt + (7^(W, y^) dW^^^ 
In detail we have 

ay{V^,yJ) = ^aJS^{Vi){l-yi)+alyix{y') (9) 
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Remember that the form of the diffusion term guarantees that the solutions to 
this equation with appropriate initial conditions stay between and 1. The 
Brownian motions W^'^ are assumed to be independent from one neuron to the 
next. 



2.4.2 Electrical synapses 

The electrical synapse transmission is rapid and stereotyped, and is mainly used 
to send simple depolarizing signals for systems requiring the fastest possible 
response. At the location of an electrical synapse, the separation between two 
neurons is very small (~ 3.5nm). This narrow gap is bridged by the gap junction 
channels^ specialized protein structures that conduct the flow of ionic current 
from the presynaptic to the postsynaptic cell (see e.g. [44j ) . 

Electrical synapses thus work by allowing ionic current to flow passively 
through the gap junction pores from one neuron to another. The usual source 
of this current is the potential difference generated locally by the action poten- 
tial. Without the need for receptors to recognize chemical messengers, signaling 
at electrical synapses is more rapid than that which occurs across chemical 
synapses, the predominant kind of junctions between neurons. The relative 
speed of electrical synapses also allows for many neurons to fire synchronously. 

We model the current for this type of synapse as 

If^'^ = Jij{t){Vi-Vj), (10) 
where Jij{t) is the maximum conductance. 



2.4.3 The maximum conductances 

As shown in equations ([6|, ^ and (10), we model the current going through 
the synapse connecting neuron j to neuron i as being proportional to the max- 
imum conductance Jij. Because the synaptic transmission through a synapse 
is affected by the nature of the environment, the maximum conductances are 
affected by dynamical random variations (we do not take into account such phe- 
nomena as plasticity). What kind of models can we consider for these random 
variations? 

The simplest idea is to assume that the maximum conductances are indepen- 

j 

dent diffusion processes with mean and standard deviation i.e. that 

depend only on the populations. The quantities J^^^ being conductances, are 
positive. We write the following equation 

where the ^ = 1, • ' ' 7 = 1? ' ' ' ^ire NP independent zero mean 

unit variance white noise processes derived from NP independent standard 
Brownian motions 5*'^(t), i.e. = ^^^it^^ which we also assume to 



11 



be independent of all the previously defined Brownian motions. This dynamics 
main advantage is its simplicity. Its main disadvantage is that if we increase 
the noise level a^^, the probability that Jij{t) becomes negative increases also: 
this would result in a negative conductance! 



One way to alleviate this problem is to modify the dynamics ( 11 ) to a slightly 
more complicated one whose solutions do not change sign, such as for instance 
the Cox-Ingersoll-Ross model ^ given by: 

dMt) = 0^,{^ - Mt))dt + ^^/j^)dB^'''{t). (12) 

Note that the righthand side only depends upon the population 7 = p{j). Let 
Jij{0) be the initial condition, it is known, [21 , that 

iE|J«(t)l = MO)'-'"' + ^(1 - «-"•-') 



This shows that if the initial condition Jij{0) is equal to the mean the 

mean of the process is constant over time and equal to Otherwise, if the 
initial condition Jij(O) is of the same sign as Jc^^, i.e. positive, then the long 
term mean is and the process is guaranteed not to touch if the condition 

2Nj6a-fJa-f ^ (^^7)^ holds Note that the long term variance is 2N3q^ • 
2.5 Putting everything together 

We are ready to write the equations of a network of Hodgkin-Huxley or Fitzhugh- 
Nagumo neurons, study their properties and their limit, if any, when the number 
of neurons becomes large. The external current for neuron i has been modelled 
as the sum of a deterministic part and a stochastic part: 



dt 



We will assume that the deterministic part is the same for all neurons in the same 
population, := /c^, and that the same is true for the variance, (7*^^ := crg^^. We 
further assume that the N Brownian motions Wl are N independent Brownian 
motions and independent of all the other Brownian motions defined in the model. 
In other words 

int) = I^(t) + a = p{i), i = l,---,N (13) 

We only cover the case of chemical synapses, and leave it to the reader to derive 
the equations in the simpler case of gap junctions. 
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2.5.1 Network of Fitzhugh-Nagumo neurons 

We assume that the parameters a^, bi and q in the equation ([5| of the adaptation 
variable w'^ of neuron i are only functions of the population a = p{i). 



Simple maximum conductances variation If we assume that the max- 



imum conductances fluctuate according to equation (11), the state of the ith 



neuron in a fully connected network of Fitzhugh-Nagumo neurons with chemical 
synapses is determined by the variables {V\w\y'^) that satisfy the following set 
of 3A^ stochatic differential equations: 



3 

P J_ 
7=1 



r 



'{t)^dt+ 



E 



i,p(i)=i "^"T 



- V,%M dt+ 



dwl 



= Ca {Vi + fla - baW]) dt 

= (a,"5„(l^/)(l - yl) - a^yj) dt + yl)dWi' 



(14) 



Scx{V^) is given by equation ([8|, by equation ([9|, and W/'^, i = 1, • • • , iV 
are N independent Brownian processes that model noise in the process of trans- 
mitter release into the synaptic clefts. 



Sign preserving maximum conductances variation If we assume that 
the maximum conductances fluctuate according to equation (12) the situation 
is slightly more complicated. In effect, the state space of the neuron i has to be 
augmented by the P maximum conductances J^^, 7 = 1, • • • , P. We obtain 



< 



^.tdWl 



dwl — ^cx iyi + <^Q! ~ bawl) dt 

dyi = (a?5«(T//)(l - yj) - a'^yj) dt + aHVlyDdWl'^ 

dJi-y{t) =9a-y{^-Ji-y{t))dt+4^^/\J~(T}\dB''''{t) J = 1,---,P 



which is a set of N{P + 3) stochastic differential equations. 



(15) 



2.5.2 Network of Hodgkin-Huxley neurons 

We provide a similar description in the case of Hodgkin-Huxley neurons. We 
assume that the functions and x e {n. m, h} that appear in equation ([2| 
only depend upon a = p{i). 



13 



Simple maximum conductances variation If we assume that the max- 
imum conductances fluctuate according to equation ( pTj ), the state of the ith 
neuron in a fuhy connected network of Hodgkin-Huxley neurons with chemical 
synapses is therefore determined by the variables (V*, n*, m*, /i*, that satisfy 
the following set of bN stochatic differential equations: 

'CdVi = (/"(t) - g-Knj{Vi - Ek) - gNamfhiiVi - E^a) - 9L{Vi - El)+) dt 
(E^=i i; T.i,,u)=f J-yi^t - V,%)yi) dt+ 

dxiit) ={p'i{y%l-Xi)-Cx{V^)xi)dt + <j^{V\Xi)dWt'' x e {n, m, /i} 
M = {a-SUVt%l - yt) - a^yt) dt + agiV^, yt)dWt " 

(16) 



Sign preserving maximum conductances variation If we assume that 



the maximum conductances fluctuate according to equation (12) we use the 
same idea as in the Fitzhugh-Nagumo case of augmenting the state space of 
each individual neuron and obtain the following set of (5 + P)N stochastic 
differential equations. 

f cdv; = (J"(t) - g-Kujiv; - Ek) - gNam^Mvi - Eno) - g-^v^ - El)+) dt 

dx.it) = - X,) - QiVi)xi) dt + a^iVl x,)dWr X G {n, m, h} 

dyi = {a^Sa{Vi){l - y\) - a^yj) dt + u^yl yt)dWr 
dJi^it) =e^^{^-Ji^{t))dt+^y^\J~(ty\dBi'-'{t) 7 = 1,---,P 

(17) 



2.5.3 Partial conclusion 



Equations (14)-(17) have a quite similar structure. They are well-posed, i.e. 
given any initial condition and any time T > they have a unique solution 
on [0, T] which is square integrable. A little bit of care has to be taken when 
choosing these initial conditions for some of the parameters, i.e. n, m and h 
which take values between and 1, and the maximum conductances when one 
wants to preserve their signs. 

In order to prepare the grounds for section |3] we explore a bit more the 
aforementioned common structure. Let us first consider the case of the simple 
maximum conductance variations for the Fitzhugh-Nagumo network. Looking 



at equation (14), we define the three-dimensional state vector of neuron i to be 
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Xi = (y/, wi, yl). Let us now define fa - ^ y.^^ , a = I,- ■ ■ ,P,\)y 



Ut,xl) 



-wi+nt) 

la - baWl) 



a?Sa{vm-yl)-a^yl 

Let us next define : R x R^^^ by 







<jyiv;,yi) 



It appears that the intrinsic dynamics of the neuron i is conveniently described 
by the equation 



dXl = Ut,Xl)dt^g^{t,Xl) 



We next define the functions 6q,^ : x ^ R^, for a, 7 = 1, • • • , P by 



0\ 








and the function (3^^ : R^ x R^ ^ R^^^ by 



Pa-f{Xl,Xl) 



It appears that the fuh dynamics of the neuron i, corresponding to equation ( 14) 
can be described compactly by 



dXi = U{t,Xi)dt + ga{t,Xi 



dWi 

dwi^y 



p 1 

Ei^ E K,{xixl)dt+ 



7=1 3,pU)=i 



7=1 



Let us now move to the case of the sign preserving variation of the maximum 
conductances, stih for Fitzhugh-Nagumo neurons. The state of each neuron is 
now P + 3-dimensional: we define XI = {V^^ wl^ yl^ Jnif)-, • • • , Jipif))- We then 
define the functions : R x R^+3 R^+^, a = 1, • • • , P, by 



a^S^{V,^){l-yl)-a2yl 
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and the functions : R x R-^+s ^ ]r(f+3) x (p+2) 



^ext 





^l{V,\yl) 




















^/\J^\ 



It appears that the intrinsic dynamics of the neuron i isolated from the other 
neurons is conveniently described by the equation 



dXi = Ut,Xi)dt + g^{t,Xi) 



dWi 

dwi'^ 



dB, 



i,l 



dBl 



i,P 



Let us finally define the functions ba^ : R^+^ x R^+^ R^+^, for a, 7 
I,-- - ,Pby 









It appears that the full dynamics of the neuron corresponding to equation ( 15 ) 
can be described compactly by 



dXi = Ut,Xi)dt^g^{t,Xi) 



dWi 

dwi;y 



dBl 



. ^7 

7=1 j,p{j)=i 



We let the reader apply the same machinery to the network of Hodgkin-Huxley 
neurons. 

We are interested in the behavior of the solutions of these equations as the 
number of neurons tends to infinity. This problem has been longstanding in 
neuroscience arousing the interest of many researchers in different domains. We 
discuss the different approaches developed in the field in the next subsection. 



2.6 Mean-field methods in computational neuroscience: a 
quick overview 

Obtaining the equations of evolution of the effective mean-field from microscopic 
dynamics is a very complex problem. Many approximate solutions have been 
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provided, mostly based on the statistical physics literature. 

Many models describing the emergent behavior arising from the interaction 
of neurons in large-scale networks have relied on continuum limits ever since 
the seminal work of Wilson and Cowan and Amari j2] |3l |82] |83| . Such mod- 
els represent the activity of the network by macroscopic variables, e.g., the 
population-averaged firing rate, which are generally assumed to be determin- 
istic. When the spatial dimension is not taken into account in the equations 
they are referred to as neural masses otherwise as neural fields. The equations 
that relate these variables are ordinary differential equations for neural masses, 
integro-differential equations for neural fields. In the second case they fall in 
a category studied in [39 or can be seen as ordinary differential equations de- 
fined on specific functional spaces [32]. Many analytical and numerical results 
have been derived from these equations and related to cortical phenomena, for 
instance for the problem of spatio-temporal pattern formation in spatially ex- 
tended models (see e.g. [20l [28l [29j [46] ) . The use of bifurcation theory has also 
proven to be quite powerful [80] [19]. Despite all its qualities, this approach 
implicitly makes the assumption that the effect of noise vanishes at the meso- 
scopic and macroscopic scales and hence that the behavior of such populations 
of neurons is deterministic. 

A different approach has been to study regimes where the activity is uncor- 
related. A number of computational studies on the integrate- and-fire neuron 
showed that under certain conditions neurons in large assemblies end up firing 
asynchronously, producing null correlations [TJ [4[ [13]. In these regimes, the 
correlations in the firing activity decrease towards zero in the limit where the 
number of neurons tends to infinity. The emergent global activity of the popu- 
lation in this limit is deterministic, and evolves according to a mean-field firing 
rate equation. However, according to the theory, these states only exist in the 
limit where the number of neurons is infinite, thereby raising the question of 
how the finiteness of the number of neurons impacts the existence and behavior 
of asynchronous states. The study of finite-size effects for asynchronous states 
is generally not reduced to the study of mean firing rates and can include higher 
order moments of firing activity [5TJ [27l [12] . In order to go beyond asynchronous 
states and take into account the stochastic nature of the firing and understand 
how this activity scales as the network size increases, different approaches have 
been developed, such as the population density method and related approaches 
[17]. Most of these approaches involve expansions in terms of the moments 
of the corresponding random variables, and the moment hierarchy needs to be 
truncated which is not a simple task that can raise a number of difficult technical 
issues (see e.g. [49]). 

However, increasingly many researchers now believe that the different intrin- 
sic or extrinsic noise sources are part of the neuronal signal, and rather than be- 
ing a pure disturbing effect related to the intrinsically noisy biological substrate 
of the neural system, they suggest that noise conveys information that can be an 
important principle of brain function [58 . At the level of a single cell, various 
studies have shown that the firing statistics are highly stochastic with probabil- 
ity distributions close to Poisson distributions [61 , and several computational 
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studies confirmed the stochastic nature of single-cehs firings [HI [57l [TH l75] . 
How does the variabihty at the single neuron level affect the dynamics of cor- 
tical networks is less well-established. Theoretically, the interaction of a large 
number of neurons that fire stochastic spike trains can naturally produce cor- 
relations in the firing activity of the population. For instance power-laws in 
the scaling of avalanche- size distributions has been studied both via models and 
experiments O |6l [48l [73] . In these regimes the randomness plays a central role. 

In order to study the effect of the stochastic nature of the firing in large 
networks, many authors strived to introduce randomness in a tractable form. 
Some of the models proposed in the area are based on the definition of a Markov 
chain governing the firing dynamics of the neurons in the network, where the 
transition probability satisfies a differential equation, called the master equation. 
Seminal works of the application of such modeling for neuroscience date back 
to the early 90s and have been recently developed by several authors [56, 15] 
[T6l [TTJ [27]. Most of these approaches are proved correct in some parameter 
regions using statistical physics tools such as path integrals and Van-Kampen 
expansions, and their analysis often involve a moment expansion and truncation. 
Using a different approach, a static mean-field study of multi population network 
activity was developed by Treves in |77] . This author did not consider external 
inputs but incorporated dynamical synaptic currents and adaptation effects. His 
analysis was completed in [1], where the authors proved, using a Fokker-Planck 
formalism, the stability of an asynchronous state in the network. Later on, 
Gerstner in [37 built a new approach to characterize the mean-field dynamics for 
the Spike Response Model, via the introduction of suitable kernels propagating 
the collective activity of a neural population in time. Another approach is based 
on the use of large deviation techniques to study large networks of neurons [33 . 
This approach is inspired by the work on spin-glass dynamics, e.g., [38 . It takes 
into account the randomness of the maximum conductances and the noise at 
various levels. The individual neuron models are rate models, hence already 
mean- field models. The mean- field equations are not rigorously derived from 
the network equations in the limit of an infinite number of neurons but they 
are shown to have a unique, non-Markov solution, i.e. with infinite memory, for 
each initial condition. 

Brunei and Hakim considered a network of integrate-and-fire neurons con- 
nected with constant maximum conductances [13^. In the case of sparse con- 
nectivity, stationarity, and in a regime where individual neurons emit spikes 
at low rate, they were able to study analytically the dynamics of the network 
and to show that it exhibits a sharp transition between a stationary regime 
and a regime of fast collective oscillations weakly synchronized. Their approach 
was based on a perturbative analysis of the Fokker-Planck equation. A simi- 
lar formalism was used in [51^ which, when complemented with self-consistency 
equations, resulted in the dynamical description of the mean-field equations of 
the network, and was extended to a multi population network. Finally, Chizhov 
and Graham [18^ have recently proposed a new method based on a population 
density approach allowing to characterize the mesoscopic behaviour of neuron 
populations in conductance-based models. 
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Let us finish tiiis very siiort and incomplete survey by mentioning the work 
of Sompolinsky and colleagues. Assuming a linear intrinsic dynamics for the 
individual neurons described by a rate model and random centered maximum 
conductances for the connections, they showed [62j |63l [22] that the system un- 
dergoes a phase transition between two different stationary regimes: a "trivial" 
regime where the system has a unique null and uncorrelated solution, and a 
"chaotic" regime in which the firing rate converges towards a non-zero value 
and correlations stabilize on a specific curve which they were able to character- 
ize. 

All these approaches have in common that they are not based on the most 
widely accepted microscopic dynamics (such as the ones represented by the 
Hodgkin-Huxley equations or some of their simplifications), and/or involve ap- 
proximations or moment closures. Our approach is distinct in that it aims 
at deriving rigorously and without approximations the mean-field equations of 
populations of neurons whose individual neurons are described by biological, if 
not correct at least plausible, representations. The price to pay is the complex- 
ity of the resulting mean-field equations. The specific study of their solutions is 
therefore a crucial step, that will be developed in forthcoming papers. 

3 Mean-field equations for conductance-based mod- 
els 

In this section, we give a general formulation of the neural networks models 
introduced in the previous section and use it in a probabilistic framework to 
address the problem of the asymptotic behavior of the networks, as the number 
of neurons N goes to infinity. In other words we derive the limit in law of N 
interacting neurons, each of which satisfying a nonlinear stochastic differential 
equation of the type described in section |2j In the remaining of this section, we 
work in a complete probability space (r2, J^, P) satisfying the usual conditions, 
and endowed with a filtration (^t)^- 

3.1 Setting of the problem 

We recall that the neurons in the network fall into P different populations. 
The populations differ through the intrinsic properties of their neurons and the 
input they receive. We assume that the number of neurons in each population 
a G {!,..., P}, denoted by A^^, increases as the network size increases, and 
moreover that the asymptotic proportion of neurons in population a is non- 
trivial, i.e. Na/N Xa ^ (O7 1) as A/" goes to infinit}|^ 

We use the notations introduced in section 12.5.31 and the reader should refer 
to this section to give a concrete meaning to the rather abstract (but required 
by the mathematics) setting that we now establish. 

■"^As we will see in the proof, most properties are valid as soon as tends to infinity as 
N goes to infinity for all o; G {1, • • • , P}, the previous assumption will allow quantifying the 
speed of convergence towards the asymptotic regime. 
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Each neuron i in population a is described by a state vector noted Xl'^ in 
and has an intrinsic dynamics governed by a drift function /q, : R x R^ R^ 
and a diffusion matrix : RxR^ ^dxm assumed uniformly locally Lipschitz 
continuous with respect to the second variable. For a neuron i in population a, 
the dynamics of the cZ-dimensional process {XI) governing the evolution of the 
membrane potential and additional variables (adaptation, ionic concentrations), 
when there is no interaction, is governed by the equation: 

dXi'"" = Ut, Xl'"") dt + g^{t, Xl'"") dWl 

We moreover assume, as it is the case for all the models described in section [2j 
that the solutions of this stochastic differential equation exist for all time. 

When included in the network, these processes interact with those of all 
the other neurons through a a set of continuous functions that only depend on 
the population a = p{i) the neuron i belongs to and the populations 7 of the 
presynaptic neurons. These functions baj{x^y) : R^ x R^ R^, are scaled 
by the coefficients so that the maximal interaction is independent of the 

size of the network (in particular, neither diverging nor vanishing as N goes to 
infinity) . 

As discussed in section [2j due to the stochastic nature of ionic currents and 
the noise effects linked with the discrete nature of charge carriers, the maxi- 
mum conductances are perturbed dynamically through the N x P independent 
Brownian motions Bl' " of dimension S that were previously introduced. The in- 
teraction between the neurons and the noise term is represented by the function 
(3aj : R^ X R^ R^x^. 

In order to introduce the stochastic current and stochastic maximum con- 
ductances, we define two independent sequences of independent m— and 6- 
dimensional Brownian motions noted (1^/)^^^ and (5^")iGN,Q;G{i - P} which are 
adapted to the filtration J^t- 

The resulting equation for the ith neuron, including the noisy interactions 
reads: 

dXl'" =Ut,Xi''')dt+f2^ E b^,{Xl'^,Xi'^)dt+g^{t,Xi''')dWi 

^ j 

El^ E P.,{Xl'\xl''')dBV. (18) 



7=1 ^ j,p{j)^l 



Note that this implies that X*'^ and X^'^ have the same law whenever p{i) = 
P{j)' 

These equations are similar to the equations studied in another context by a 
number of mathematicians, among which McKean, Tanaka and Sznitman (see 
the Introduction), in that they involve a very large number of particles (here 
particles are neurons) in interaction. Motivated by the study of McKean- Vlasov 



equations, these authors studied special cases of equations (18). This theory 



20 



referred to as the kinetic theory, is chiefly interested in the study of the ther- 
modynamics questions. They show the property that in the hmit where the 
number of particles tends to inflnity, provided that the initial state of each 
particle is drawn independently from the same law, each particle behaves in- 
dependently and has the same law, which is given by an implicit stochastic 
equation. They also evaluate the fluctuations around this limit, under diverse 
conditions |^ EH [TTl [65j [66] [53l [52] . Some extensions to biological problems 
where the drift term is not globally Lipschitz but satisfles the monotone growth 



condition (19) were studied in [7]. This is the approach we undertake here. 
We now spell out our assumptions on the function appearing in the model: 

(HI). Locally Lipschitz dynamics: for all ol G {1,...,P}, the functions /q, 
and Qa, are uniformly locally Lipschitz-continuous with respect to the sec- 
ond variable. In detail, for all > there exists Ku > and independent 
of t e [0, T] such that for ah x, y e Bfj, the bah of of radius U\ 

\\fa{t, X) - fa{t, y)\\ + \\ga{t, x) - ga{t, ^) || < i^f/ ||x - ^|| a = 1, • • • , P, 

(H2). Lipschitz interactions: for all 0^,7 G {1, . . . ,P}, the functions ba-f and 
I3a-f are Lipschitz-continuous with Lipschitz constant L: for all (x, y) and 
(x^ y^) in x R^ we have: 

\K^{x,y)-ba^{x\y')\^\\(3a^{x,y)-(3a^{x\y'^^^ 

(H3). Linear growth of the interactions: There exists a K > such that: 
max(||6„^(x,z)f , \\p^^{x,z)f) < k{l + \\xf) 

(H4). Monotone growth of the dynamics: We assume that fa and satisfy 
the following monotonous condition for all a = 1, • • • , P: 

x^Ut.x) + l\\ga{t^x)f <K{1^ \\xf) (19) 

3.2 Convergence of the network equations to the mean- 
field equations and properties of those equations 

We now show that the same type of phenomena that were predicted for systems 
of interacting particles happen in networks of neurons. In detail we prove that 
in the limit of large populations the network displays the property of propaga- 
tion of chaos. This means that any finite number of diffusion processes become 
independent, and all neurons belonging to a given population a have asymptot- 
ically the same probability distribution, which is the solution of the following 
mean-field equation: 



dX^ = Uit, Xn dt + J2 ^zKyiXr, ZD\ dt + gUt, Xn dW^ 

p 

Y,^^[Pa,{K.mdBr a = l,---,P (20) 



7=1 

p 



7=1 
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where Z is a process independent of X that has the same law, and denotes 
the expectation under the law of Z. In other words, the mean-field equation 
can be written, denoting by m]{dz) the law of (hence also of X^): 

7=1 ^-^R" 



z)dnq{z)\dBr (21) 



In these equations, ^ for a = 1 • • • P, are independent standard m-dimensional 



Brownian motions. Let us point out the fact that the righthand side of (20) and 



(21 ) depends on the law of the solution, this fact is sometimes referred to as "the 
process X is attracted by its own law" . This equation is also classically written 
as the McKean-Vlasov-Fokker-Planck equation on the probability distribution 
p of the solution. This equation which we use in section [4] can be easily derived 
from equation (20). Let Pa{t^ z)^ z = (^i, • • • , Zd)^ be the probability density at 



time t of the solution X^ to (20) (this is equivalent to dmf{z) = Pa{t^ z)dz), we 
have: 



dtPa{t,z) 



div. 



^ r 



iz,y)Pj it,y) dy {t,z) 



d 



2 ^ dzidzj 



{Dfj{z)p^{t,z)) 



a 



1,--,P, (22) 



where the d x d matrix is given by 

p 



D^{z) = Y,^z[Pa^{z,Z)]^l[p^^{z,Z)]. 



7=1 



with 



y{z,Z)]= J Pa-f{z,y)p^{t,y)dy. 



The P equations (22) yield the probability densities of the solutions of 
the mean-field equations (20). Because of the propagation of chaos result, the 



are statistically independent, but their probability functions are clearly 
functionally dependent. 

We now spend some time on notations in order to obtain a somewhat more 



compact form of equation ( 20 ) . We define Xf to be the d P-dimensional process 
Xf = (Xf ; a = 1 • • • P). We similarly define f^g^b and P the concatenation of 
functions ga, ba,f3 and In details f{t,X) = {fa{t,X^) ; a = 1 • • - P), 

b{X,Y) = iE^^iba^iX-.Y^) ; a = 1---P)), and W = {W^ ; a = 1---P). 
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The term of noisy synaptic interactions requires a more careful treatment. We 



define P = {13 



a, 7 



dx5 



PxP 



and B = (5° 



and the product of an element M G (R' 



dx6 



PxP 



element X G (R^)^^^ the element of (R^^) 

(M0X), = ^M,^X^ 



d\P. 



a,7 = 
and an 



We obtain the equivalent compact mean-field equation: 
dXt= (^f{t,Xt)^]Ez[b{Xt,Zt)]^dt^g{t,Xt)dWt^lEz (23) 



The equations (20) and (22) are implicit equations on the law of Xt 



We now state the main theoretical results of the paper as two theorems. 



The first theorem is about the well-posedness of the mean-field equation (20). 



The second is about the convergence of the solutions of the network equations 
to those of the mean-field equations. Since the proof of the second theorem 
involves similar ideas to those used in the proof of the first, it is given in the 
appendix [A) 



Theorem 2. Under the assumptions ] (Hl)\ to \( H4 )\ there exists a unique solution 
to the mean-field equation (po]) on [0,T] for any T > 0. 



Let us denote by A4{C) the set of probability distributions on C the set con- 
tinuous functions [0,T] (R^) , and M'^{C) the space of square-integrable 
processes. Let also (W" ; a = 1 • • • P) (respectively (P^^, a,7 = 1---P)) 
be a family of P (respectively P^) independent m- (respectively (5)-dimensional 
adapted standard Brownian motions on (Q^J^^P). Let us also note Xq G 
A1(R^)^ the (random) initial condition of the mean-field equation. We in- 
troduce the map ^ acting on stochastic processes and defined by: 



^ : <^ 



(M{C) 
X 



^(r, = {rf,a 



1---P})t with 



X^ 



Jo {fa{s, xf) + ^zK^{xf, Z2)]) ds^ 

f^g^{s,X-)dW- 



I,--- ,p 



We have introduced in the previous formula the process Zt with the same law 
as and independent of Xf. There is a trivial identification between the solutions 



of the mean-field equation (20) and the fixed points of the map $: any fixed 
point of ^ provides a solution for equation (20) and conversely any solution of 



equation (20) is a fixed point of ^. 



The following lemma is useful to prove the theorem. 

Lemma 3. Let Xq G ]L^((R^)^) a square-integrable random variable. Let X be 
a solution of the mean-field equation (20) with initial condition Xq. Under the 
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assumption (19), there exists a constant C{T) > depending on the parameters 
of the system and on the horizon T, such that: 

nUt\?\<C{T), \/te[Q,T] 

Proof. Using Ito formula for we have: 

llX^f = ||Xo||2 + 2^* (xJ/(s,X,) + i||5(s,X,)f + XjEz[6(X„Z,)] + 

i||Ez[/3(X„Z,)]f))ds + iV, 
where Nt is a stochastic integral, hence with a null expectation, ]E[A^^] = 0. 



This expression involves the term x^b{x,z). Because of assumption (H3) 
we clearly have: 



\x^b{x,z)\ < \\x\\ \\b{x,z)\\ < \\x\\^K{l + \\x\\^) < ^/R{l+\\xf) 
It al so inv olves the term x^ f{t^x) + which, because of ass ump- 



tion (H4), is upperbounded by K{1 + H^^lP). Finally, assumption (H3) 
allows us to upperbound the term ^\\]Ez[P{Xs, Zs)]f by f (1 + \\Xsf). 
Finally we obtain 



agam 



E[l + \\Xtf] = E[l + \\Xof] +2{K+^ + ^) J^n^+ \\Xsf] ds. 

Using Gronwall's inequality, we deduce the boundedness of the solutions of 
the mean- field equations. □ 

This lemma puts us in a position to prove the existence and uniqueness 
theorem: 

Proof We start by showing the existence of solutions, and then prove the 
uniqueness property. We recall that by application of lemma |3| the solutions 
will all have bounded second order moment. 
Existence: 

Let = {X^ = {X^""^ a = 1---P}) e M{C) be a given stochastic process, 
and define the sequence of probability distributions {X^)k>o on M{C) defined 
by induction by = ^{X^). Define also a sequence of processes Z^, /c > 0, 

independent of the sequence of processes X^ and having the same law. We note 
this as "X and Z i.i.d." below. We stop the processes at the time the first 
hitting time of the norm of X^ to the constant value U. For convenience, we will 
make an abuse of notation in the proof and denote Xj^ = X^^^^ • This implies 

that X^ belongs to for ah times t G [0,T]. 
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Using the notations introduced for equation ( 23 ) we decompose the difference 
X^^^ - as fohows: 



X^^' - X^ = 1^ [f{s,X^) - f{s,Xt')) 



ds 



At 



[\z\b{XlZ^)-b{Xt\Zt'] 
Jo L 



ds 



Bt 



£ (g{s,X^,)-g{s,X^,-'))dWs 



f 

Jo 



Ct 



/3(x^z,^)-/3(x^^z^l) 



(ddBs 



Dt 



and find an upperbound for 



k ._ 



E 



sup,<^ \\X^^^ - X^W^ by finding up- 



perbounds for the corresponding norms of the four terms A^, Bt, Ct and Dt. 
Applying the discrete Cauchy- Schwartz inequahty we have: 



-fe+l _ Y-/c||2 ^ ^ 



{\\At 



\\BA\^ 



\Ct\ 



and treat each term separately. The upperbounds for the first two terms are 
obtained using the Cauchy- Schwartz inequality, those of the last two terms using 
the Burkholder-Davis-Gundy martingale moment inequality. 

The term At is easily controlled using the Cauchy- Schwarz inequality. In 
detail 



(Cauchy- Schwartz) < T 



r {f{u,Xt)-f{u,Xt')) du 
Jo 

^ r\\f{u,Xt)-f{u,Xt')\f du 
Jo 

(assumption pT)] ) < K^T [ \\X^ - X^-^f du 

Jo 

Taking the sup of both sides of the last inequality we obtain 
sup||A,||2 < kIjT f -X^^ll' ds < KlT f suv\\xl-Xl-^f ds, 

s<t Jo Jo u<s 



from which follows the fact that 



E 



sup 11^5 

s<t 



< KjjT 



E 



sup ||X^ 



X, 



fc-ii 



ds 
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The term Bt is controlled using the Cauchy- Schwartz inequality, the uniform 
continuity of the function b and the fact that the processes X and Z are inde- 
pendent with the same law. In detail 

llS.f = riEz[b{Xt,Zt)-b{Xt\Zt-')] du 
Jo 

(Cauchy-Schwartz) <T [ ||Ez [b{X^, Z^) - b{X^-\ Z^'^)] f du 
Jo 



(Cauchy-Schwartz) < T 



b{xlzt)-b{x 



k—l yk — \\ 

U 1 ) 



du 



(assumption [(H2)} < 2TLM Ez \\X^ - X^ 



Z' k /yk 



7k-\ I 



du 



{X and Z i.i.d.) < 2TL^ 



(II 



Taking the sup of both sides of the last inequality we obtain 



sup 115,11^ <2Ti:^ / [Ws-X, 

s<t 



1> 



k-1 1 



•E 



X^-X, 



k-l I 



ds 



<2TL' [ fsup||X^-X^ 
^0 V 



k-l I 



- supE 



-k-l I 



Xk \rk 



<2TL^ [ (sup||X„^-X^l| 

Jo \u<s 



-E 



sup||X^-X, 



k-l\ 
u ^^u I 



ds 
ds, 



from which follows the fact that 



E 



sup \\Bs 

.s<t 



<4TL^ / E 



/ 

JO 



sup \\Xl - X, 

u<s 



k-l I 



ds 



The term Ct is controlled using the fact that it is a martingale and applying 
the Burkholder-Davis-Gundy martingale moment inequality: 









E 


sup||C,f 


<4 / E 




_S<t 


Jo ^ 



\g{s,X^)-9{s,Xt')f 



ds 



(assumption [(HI)] ) < AK^ [ E - X^-^\ 

Jo ^ 



ds 



<4K^ I E 



sup\\X^-Xt'\f 



u<s 



ds 



The term Dt is also controlled using the fact that it is a martingale and applying 
the Burkholder-Davis-Gundy martingale moment inequality: 
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E 



sup II A f 



,s<t 



< 4 / E 



(Cauchy-Schwartz) < 4 / E 



(assumption |(H2)| ) < 8L'^ [ E 

Jo 

{X and Z i.i.d.) < WL^ [ E 
Jo 

< 16L^ [ E 
Jo 

Putting all this together, we get: 



Ez 
Ez 



x^,z^)-P{xr\z, 



ds 
ds 



ds 



ds 



sup||X„«-X„« 



k-1 1 



E 



sup||Xr^-X,^||^l <4(T + 4)(i^^+4L^) / E[sup - 



u<s 



(24) 



From the relation < K" M^'^ ds with K" = 4(T + 4)(i^2 +4^^), we 
get by an immediate recursion: 



Jo JO Jo 



fcT 



< 



(25) 



and is finite because the processes are bounded. Bienayme-Tchebychev 
inequality and (25) now give: 



< 4^^ — tt^MS. 



ans this upper bound is the term of a convergent series. From the Borel-Cantelli 
lemma stems that for almost any u e Q there exists a positive integer ko{uj) {uj 
denotes an element of the probability space Q) such that 



sup ||X, 

S<t 



/c+1 



and hence 



22(/c+l) 
1 



sup|ixr^-x,^ii<2,^^ 



V A: > ko{uj), 



V/c > ko{oj). 



It follows that with probability 1 the partial sums: 

n 

X° + ^(X^i - X^) = X- 



k=0 
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are uniformly (in t G [0,T]) convergent. Denote by Xf the thus defined hmit. 
It is clearly continuous and Tt adapted. On the other hand, the inequality (25) 
shows that for every fixed the sequence {X^}n>i is a Cauchy sequence in L . 
Lemma [3] shows that X e M^{C). 

It is easy to show using routine methods that X indeed satisfies the equa- 
tion ( [2Q1 ). 

To complete the proof we use a standard truncation property. This method 
replaces the function / by the truncated function: 



fu{t,x) 



f(t,x) \\x\\<U 
f{t,Ux/\\x\\) \\x\\>U 



and similarly for g. The functions fu and gu are globally Lipchitz-continuous, 
hence the previous proof shows that there exists a unique solution Xu to equa- 
tions (20) associated with the truncated functions. This solution satisfies the 
equation 

Xu(t) = Xo+ f {fu{t,Xu{s)) ^1^2[KXu{s),Z,)]) ds^ f gu(t,Xu{s))dWs 
Jo Jo 

+ [ lEz [l3{Xu{s), Z,)] QdB, te [0, T] (26) 
Jo 

Let us now define the stopping time: 

Tu=ini{tG[0,T],\\Xum>U}. 

It is easy to show that 

Xu{t)=Xu'{t) if 0<t<ru,U'>U, (27) 

implying that the sequence of stopping times tu is increasing. Using lemma [S] 
which implies that the solution to (20) is almost surely bounded, for almost 
all 00 e ft^ there exists Uo{uj) such that ru = T for all U > Uq. Now define 
X{t) = Xuoit), t ^ [0, T]. Because of ^ we have X{tATu) = Xu{t/\ru) and 
it follows from (26) that 

X{tATu) = Xo + (^fu{s,Xs)+^z[b{Xs,Zs)])ds + gu{s,Xs)dWs 

ptAru 

+ / ]Ez[P{Xu{s),Zs)]&dBs 
Jo 



ptAru 

+ / ]E2[P{Xu{s),Zs)]QdBs 
Jo 



and letting U ^ oo we have shown the existence of solution to equation (20) 
which by lemma [3] is square integrable. 
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Uniqueness: 

Assume that X and Y are two solutions of the mean-field equations (20). From 
lemmajsj we know that both solutions are in A^^(C). Moreover, using the bound 
(24), we directly obtain the inequality: 



E 



sup||x,-r,| 



Jo 



E 



sup \\Xu - Yu\ 

u<s 



ds 



which by Gronwall's theorem directly implies 



E 



sup \\Xs - Ys 

s<t 



which ends the proof. 



□ 



We have proved the well-posedness of the mean-field equations. It remains 
to show that the solutions to the network equations converge to the solutions 
of the mean-field equations. This is what is achieved in the next theorem. 

Theorem 4. Under the assumptions \(H1)\ to \(H4)\ the following holds true: 

• Convergenc^ For each neuron i of population a, the law of the multi- 
dimensional process X*'^ converges towards the law of the solution of the 
mean-field equation related to population a, namely X" . 

• Propagation of chaos: For any k G N*^ and any k-uplet (zi, . . . , ij^), the 
law of the process (X^^'^, . . . ,X^*'^'^,t < T) converges toward^m^ 
... m^^*'^^ i.e. the asymptotic processes have the law of the solution of 
the mean-field equations and are all independent. 

This theorem has important implications in neuroscience that we discuss in 
section |5] Its proof is given in Appendix |Aj 



4 Numerical simulations 

At this point, we have provided a compact description of the activity of the 
network when the number of neurons tends to infinity. However, the struc- 
ture of the solutions of these equations is complicated to understand from the 



implicit mean-field equations (20) and of their variants (such as the McKean- 
Vlasov-Fokker-Planck equations (|22|)). In this section, we present some classical 
ways to numerically approximate the solutions to these equations and give some 
indications about the rate of convergence and the accuracy of the simulation. 
These numerical schemes allow us to compute and visualise the solutions. We 
then compare the results of the two schemes for a network of Fitzhugh-Nagumo 
neurons belonging to a single population and show their good agreement. 

^The type of convergence is specified in the proof given in Appendix [a| : 
^The notation mf was introduced right after equation \20\ . 
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The main difficulty one faces when developing numerical schemes for equa- 



tions (20) and (22) is that they are non local. By this we mean that, in the case 



of the McKean-Vlasov equations, they contain the expectation of a certain func- 



tion under the law of the solution to the equations, see (20). In the case of the 



corresponding Fokker-Planck equation, it contains integrals of the probability 



density functions which is solution to the equation, see (22). 



4.1 Numerical simulations of the McKean-Vlasov equa- 
tions 

The fact that McKean-Vlasov equations involve an expectation of a certain 
function under the law of the solution of the equation makes them particularly 
hard to simulate directly. One is often reduced to use Monte-Carlo simula- 
tions to compute this expectation, which amounts to simulating the solution of 
the network equations themselves (see [68]). This is the method we used. In 
its simplest fashion, it consists of a Monte-Carlo simulation where one solves 



numerically the N network equations (18) with the classical Euler-Maruyama 



method a number of times with different initial conditions, and averages the 
trajectories of the solutions over the number of simulations. 

In detail let At > and N G N*. The discrete-time dynamics implemented 
in the stochastic numerical simulations consists in simulating M times a P- 
population discrete-time process < T/At^i = 1---7V), solution of the 

recursion, for i in population a: 



7=1 ^ J = l,p(i)=7 

+ v^{^«(t,xr)C';i + EAr E /^«7(^r,^r)-c+i} (28) 

7=1 ^ J = l,p(j)=7 

where ^^'^ and Cn^'^ are independent d- and (^-dimensional standard normal 
random variables. The initial conditions i = ,A/' are drawn in- 

dependently from the same law within each population for each Monte-Carlo 
simulation r = 1, • • • , M. One then chooses one neuron in each population 
a = 1, • • • , P. If the size TV of the population is large enough, theorem [4] states 
that the law, noted Pa{t^X)^ of X*" should be close to that of the solution X" 
of the mean-field equations, for a = 1, • • • , P. Hence, in effect, simulating the 
network is a good approximation, see below, of the simulation of the mean- 
field or McKean-Vlasov equations [8, 68 . An approximation of Pa{t,X) can be 
obtained from the Monte-Carlo simulations by quantizing the phase space and 
incrementing the count of each bin whenever the trajectory of the neuron 
at time t falls into that particular bin. The resulting histogram can then be 



compared to the solution of the McKean-Vlasov- Fokker-Planck equation (22) 
corresponding to population a whose numerical solution is described next. 
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The mean square error between the solution of the numerical recursion ( 28 ) 
and the solution of the mean-field equations (20) X^/^t order 0(VAt + 



1/vA^), the first term being related to the error made by approximating the 
solution of the network of size TV, X^'^^ by an Euler-Maruyama method, and 
the second term to the convergence of X^^^ towards the mean-field equation 
^nAt when considering globally Lipschitz-continuous dynamics (see proof of 
theorem |4] in appendix [A]) . 

4.2 Numerical simulations of the McKean-Vlasov-Fokker- 
Planck equation 



For solving the McKean-Vlasov-Fokker-Planck equation ([22| we have used the 
method of lines [59l |60] . Its basic idea is to discretize the phase space and 
to keep the time continuous. In this way the values Pa{t,X), a = !,••• ,P 
of the probability density function of population a at each sample point X of 
the phase space are the solutions of P ordinary differential equations where 
the independent variable is the time. Each sample point in the phase space 
generates P ODEs, resulting in a system of coupled ODEs. The solutions to this 



system yield the values of the probability density functions Pa solution of (22) 
at the sample points. The computation of the integral terms that appear in the 
McKean-Vlasov-Fokker-Planck equation is achieved through a recursive scheme, 
the Newton-Cotes method of order 6 [79 . The dimensionality of the space being 
large and numerical errors increasing with the dimensionality of the integrand, 
such precise integration schemes are necessary. For an arbitrary real function / 
to be integrated between the values xi and X2, this numerical scheme reads: 

/ f{x)dx ^ Ax y [19/(xi + {5i - 5) Ax) + 75/(xi + {5i - 4) Ax)+ 
Jx, 288 ^ 

+50/(xi+(5i-3)Ax)+50/(xi+(5z-2)Ax)+75/(xi+(5i-l)Ax)+19/(xi+5iAx)] 

where Ax is the integration step and M = {x2 — Xi)/Ax is chosen to be an 
integer multiple of 5. 

The discretization of the derivatives with respect to the phase space param- 
eters is done through the following fourth-order central difference scheme: 

df{x) _ f{x - 2Ax) - Sf{x - Ax) + Sf{x + Ax) - f{x + 2Ax) 
~dx 12Ax ' 

for first-order derivatives, and 

d'^f{x) _ -f{x - 2Ax) + 16f{x - Ax) - 30f{x) + 16f{x + Ax) - f{x + 2Ax) 
dx^ ^ 12Ax^ ' 

for the second-order derivatives, see [54 . 

Finally we have used a Runge-Kutta method of order 2 (RK2) for the nu- 
merical integration of the resulting system of ODEs. 
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4.3 Comparison between the solutions to the network and 
the mean-field equations 

As advertised, we illustrate these ideas with the example of a network of 100 
Fitzhugh-Nagumo neurons belonging to one, excitatory, population. We choose 



a finite volume outside of which we assume that the p.d.f., solution of (22), is 
zero. We then discretize this volume with nynw^v points defined by 

nv =^ {Vmax - Vmin)/^V 
{Umax ymin) / 

where Vmin, Vmax, Wmin, Wmax, Vmin, Vmax define the volumc in which we solve 
the McKean-Vlasov-Fokker-Planck equation and estimate the histogram defined 



in section 4.1 while AV, Ak;, A?/ are the quantization steps in each dimension 
of the phase space. 

In general, the total number of coupled ODEs is the product Pnynwriy (in 
our case we chose P = 1). This can become fairly large if we increase the 
precision of the phase space discretization. Moreover, increasing the precision 
of the simulation in the phase space requires, in order to ensure the numerical 
stability of the method of lines, to decrease the time step At used in the RK2 
scheme. This can strongly impact the efficiency of the numerical method, see 
section 14.41 

In the simulations shown in the lefthand parts of figures |4] and [5] we have 
used one population of 100 excitatory FitzHugh-Nagumo neurons connected 
with chemical synapses. We performed 10, 000 Monte Carlo simulations of the 



network equations ( 14 ) with the Euler-Maruyama method in order to approxi- 
mate the probability density. The model for the time variation of the synaptic 
weights is the simple model. The p.d.f. p(0, y) of the initial condition is 
Gaussian and reads 

(V_Vq)2 (^-^o)2 (y-^o)2 

p(0,F,^,y) = ----^ e ^^-0 ^< '-lo (29) 

The parameters are given in the first column of table [l] In this table, the 
parameter tfin is the time at which we stop the computation of the trajecto- 
ries in the case of the network equations, and the computation of the solution 
of McKean-Vlasov-Fokker-Planck equation in the case of the mean-field equa- 
tions. The sequence [0.5, 1.5, 1.8, 3.0] indicates that we compute the solutions 
at those four time instants corresponding to the four rows of figure [5] The means 
take two different values corresponding to the two numerical experiments de- 
scribed below. The phase space has been quantized with the parameters shown 
in the second column of the same table to solve the Fokker-Planck equation. 
This quantization has also been used to build the histograms that represent 
the marginal probability densities with respect to the pair (F, y) of coordinates 
of the state vector of a particular neuron. These histograms have then been 
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interpolated to build the surfaces shown in the lefthand side of figures |4]and[5l 
The parameters of the Fitzhugh-Nagumo model are the same for each neuron 
of the population: they are shown in the third column of table [T] Two values of 
the external current are given, corresponding to the two numerical experiments 
whose results are shown below. 

The parameters for the noisy model of maximum conductances of equa- 
tion ( 11 ) are shown in the fourth column of the table. For these values of J and 
(jj the probability that the maximum conductances change sign is very small. 
Finally, the parameters of the chemical synapses are shown in the sixth column. 
The parameters F and A are those of the x function ([3|. The solutions are 
computed over an interval of t fin = 3 time units for figure [4] and for the four 
values shown in brackets for figure [5] with a time sampling of At = 0.1. 



The marginals estimated from the trajectories of the network solutions are 
then compared with those obtained from the numerical solution of the McKean- 
Vlasov-Fokker-Planck equation (see figures |4] and [s] right), using the method of 
lines explained above, and starting from the same initial conditions (29) as the 
neural network. 



Initial Condition 


Phase space 


FitzHugh- 
Nagumo 


Synaptic 
Weights 


Synapse 


tfin — 3.0, 


^rain 3 


a = 0.7 


J = 1 


Vrev 1- 


[0.5,1.5,1.8,3.0] 










At = 0.1 


^max — 3 


6 = 0.8 


a J = 0.2 


ar = 1 


Fo = -1.472, 0.0 


AV = 0.1 


c = 0.08 




ad = l 


avo = 0.2 




/ = 0, 0.7 




T — 1 

max 


^0 = -0.965, 0.5 


^max 2 


O'ext 




A = 0.2 


= 0-2 


A^ = 0.1 






Ft = 2 


^0 = 0.250, 0.3 


l/min ~ 






r = o.i 


ay, = 0.05 


Umax — 1 

Ay 0.06 






A = 0.5 



Table 1: Parameters used in the simulations of the neural network and for 
solving the McKean-Vlasov-Fokker-Planck equation whose results are shown in 
figures [4]and|5l see text. 



In a first series of experiments we choose the external input current value to 
be / = (see third column of table [T]). This corresponds to a stable fixed point 
for the isolated Fitzhugh-Nagumo neuron. 

The initial conditions (see first column of table [T]) are chosen in such a way that 
on average, the initial point is very close to this fixed point. We therefore expect 
that the solutions of the neural network and the McKean-Vlasov-Fokker-Planck 
equation will not change much over time and will concentrate their mass around 
this fixed point. This is what is observed in figure |4j where the simulation of 
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the neural network (lefthand side) is in excellent agreement with the results of 
the simulation of the McKean-Vlasov- Fokker-Planck equation (righthand side). 



Fokker-Planck Equation 




Figure 4: Joint probability distribution of (V, ?/), (Left): Monte-Carlo simulation 
of the McKean-Vlasov equations (20) and (Right): Simulation of the McKean- 
Vlasov- Fokker-Planck equation (22). Parameters of table [l] / = 0. These 
parameters corresponds to a stable equilibrium. Initial conditions (first column 
of table [T]) are concentrated around this point. The two distributions are very 
similar and centered around the equilibrium point, see text. 



We then changed the value of the external current to / = 0.7 (this value 
corresponds to the existence of a stable limit cycle for the isolated Fitzhugh- 
Nagumo neuron), keeping the other parameters values the same except for the 
initial conditions which now take the values Vq = 0, wq = 0.5, = 0.3. We 
have also chosen four values of tfin^ i.e. 0.5, 1.5, 1.8 and 3.0, instead of only 
one in the previous numerical experiment. They correspond to four snapshots 
of the network dynamics. We therefore expect that the solutions of the neu- 
ral network and the McKean-Vlasov-Fokker-Planck equation will concentrate 
their mass around the limit cycle. This is what is observed in figure [5j where 
the simulation of the neural network (lefthand side) is in excellent agreement 
with the results of the simulation of the McKean-Vlasov-Fokker-Planck equa- 
tion (righthand side). Note that the densities display two peaks. These two 
peaks correspond to the fact that, depending upon the position of the initial 
condition with respect to the nullclines of the Fitzhugh-Nagumo equations, the 
points in the phase space follow two different classes of trajectories, as shown in 
figure [6] The two peaks then rotate along the limit cycle, see also section [44l 



4.4 Numerical simulations with GPUs 

Unfortunately the algorithm for solving the McKean-Vlasov-Fokker-Planck equa- 
tion described in the previous section is computationally very expensive. In fact, 
when the number of points in the discretized grid of the (V^w^y) phase space is 
big, i.e. when the discretization steps AV^ Aw and Ay are small, we also need 
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Fokker-Planck Equation 



Neural Network 




Potential V 2 3 0.0 

Fokker-Planck Equation 



Neural Network 



Potential V 

Fokker-Planck Equation 



Conductance Y 




Conductance Y 



Neural Network 




Figure 5: Joint probability distribution of (F, y) computed with the Monte- Carlo 
algorithm for the network equations (20) (left) compared with the solution of 
the McKean-Vlasov- Fokker-Planck equation ( [22| ) (right), sampled at for four 
times tfin. Parameters given in table [l] with a current / = 0.7 corresponding to 
a stable limit cycle. Initial conditions (first column of table [T]) are concentrated 
inside this limit cycle. The two distributions are similar and centered around 
the limit cycle with two peaks, see text. 
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Potential V 



Limit Cycle in the W-Y piane 




Adaptation W 



Figure 6: Projection of 100 trajectories in the (F, w) (top left), (V, y) (top right) 
and {w, y) (bottom) planes. The limit cycle is especially visible in the (F, w) 
projection (red curves). The initial conditions split the trajectories into two 
classes corresponding to the two peaks shown in figure [5] Parameters similar to 
those of figure [5] 



36 



to keep At small enough in order to guarantee the stability of the algorithm. 
This implies that the number of equations that must be solved has to be large 
and moreover that they must be solved with a small time step if we want to 
reduce the numerical error. This will inevitably slow down the simulations. We 
have dealt with this problem by using a more powerful hardware, Graphical 
Processing Units (GPUs). CPUs are low cost highly parallel devices which are 
currently not only used for image or video processing but also for numerical 
simulations in physics or other scientific areas. They are specially designed 
for running a huge amount of low cost threads featuring a high speed context 
change between them. All these capabilities make them suitable for solving the 
McKean-Vlasov-Fokker-Planck equations. 

We have improved the Runge-Kutta scheme of order 2 used for the simula- 
tions shown in section [4j and adopted a Runge-Kutta scheme of order 4. 

One of the purposes of the numerical study is to get a feeling for how the 
different parameters, in particular those related to the sources of noise, influence 
the solutions of the McKean-Vlasov-Fokker-Planck equation. This is meant to 
prepare the ground for the study of the bifurcation of these solutions with re- 
spect to these parameters, as was done in [76] in a different context. For this 
preliminary study we varied the input current / and the parameter aext con- 



trolling the intensity of the noise on the membrane potential in equations (14) 
The McKean-Vlasov-Fokker-Planck equation writes in this cas^ 



d 

—p{t,V,w,y) 

-I 



3 



V-Y- W+I+J{V-Vr, 



,) [ y'p{t,V',w',y')dV'dw'dy']p{t,V,w,y)} 



— [c(V + a - bW)p{t, V, w, y)] - —{[arS{V){l - y) - ady]p{t, V, w, y)}+ 



1 32 , 2 



ext 



■ajiV-Vrevf I y'p{t,V\w\y')dV'dw'dy']p{t,V,w,y)] 



2dV^ 

\<^l^P{t^y^w^y) + \^{WS{V){l-y)+aay]x^{y)p{t,V,w,y)} (30) 

The simulations were run with the x function ([3|, the initial condition de- 
scribed by ( 29 ) and with the parameters shown in table[2j Three snapshots of the 



solution are shown in figures [t] (corresponding to the values / = 0.4, (Text = 0-27 
of the external input current and of the standard deviation of the noise on the 
membrane potential) and [s] (corresponding to the values / = 0.7, (Text = 0.45). 
In the figures the left column corresponds to the values of the marginal p(t, V, w)^ 
the right column to the values of the marginal pit.V^y). Both are necessary to 
get an idea of the shape of the full distribution pit^V^w^y). The first row of 



^We have included a small noise (controlled by the parameter (Tw) on the adaptation 
variable w. This does not change the previous analysis, in particular proposition^ but makes 
the McKean-Vlasov-Fokker-Planck equation well-posed in a cube of the state space with 
boundary value, see e.g., t31j . 
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Figure 7: Marginals with respect to the V and w variables (Left) and to the 
V and y variables (Right) of the solution of the McKean-Vlasov-Fokker-Planck 
equation. The first row shows the initial conditions, the second the marginals at 
time 0.05 and the third the stationary (large time) solutions. The input current 
/ is equal to 0.4 and (Text = 0.27. These are screenshots at different times of 
movies available as supplementary material. The corresponding movies are in 
the files Vw-n0271-in04.avi and Vy-n0271-in04.avi. 
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Figure 8: Marginals with respect to the V and w variables (Left) and to the V 
and y variables (Right) of the solution of the McKean-Vlasov-Fokker-Planck 
equation. The first row shows the marginals at time 0.05, the second the 
marginals at time 0.10 and the third the stationary (large time) solutions. The 
input current / is equal to 0.7 and (Text = 0.45. These are screenshots at dif- 
ferent times of movies available as supplementary material. The corresponding 
movies are in the files Vw-n0451-in07.avi and Vy-n0451-in07.avi. 
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figure [7| shows the initial conditions. They are the same for the results shown in 
figure [8) The second and third rows of figure [7| show the time instants t = 0.05 
and at convergence t = 0.33 (the time units differ from those of the previous 
section but it is irrelevant to this discussion) . The three rows of figure [s] show 
the time instants t = 0.05, t = 0.10 and t = 0.33. In both cases the solution 
appears to converge to a stationary distributions whose mass is distributed over 
a "blurred" version of the limit cycle of the isolated neuron. The "blurriness" 
increases with the variance of the noise. The four movies for these two cases are 
available as supplementary material. 
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a J = 0.01 


Vo = 0.0 


AV = 0.027 


c = 0.08 




avo = 0.2 
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wq = —0.5 


'^max — 3 


(Je^t = 0.27,0.45, 




crn,o = 0.2 


Aw = 0.02 


cr^ 0.0007 




Vo = 0.3 


ymin — 






ay„ = 0.05 


Umax — 1 








= 0.003 







Table 2: Parameters used in the simulations of the McKean-Vlasov-Fokker- 
Planck equation with GPUs shown in figures [7] and [Sj and in the supplementary 
material. 



5 Discussion and conclusion 

In this article, we addressed the problem of the limit in law of networks of 
biologically-inspired neurons as the number of neuron tends to infinity. We em- 
phasized the necessity of dealing with biologically-inspired models and discussed 
at length the type of models relevant to this study. We chose to address the case 
conductance-based network models that are relevant description of the neuronal 
activity. Mathematical results on the analysis of these diffusion processes in in- 
teraction resulted to replace a set of NP (i-dimensional coupled equations (the 
network equations) in the limit of large A^s by P d-dimensional mean-field equa- 
tions describing the global behavior of the network. However, the price to pay 
for this reduction was the fact that the resulting mean-field equations are non- 
standard stochastic differential equations, similar to McKean-Vlasov equations. 
These can either be expressed as implicit equations on the law of the solution, 
or, in terms of probability density function through the McKean-Vlasov-Fokker- 
Planck equations as a nonlinear, nonlocal partial differential equation. These 
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equations are in general hard to study theoretically. 

Beside the fact that we explicitely model real spiking neurons, the mathe- 
matical part of our work differs from that of previous authors such as McKean, 
Tanaka and Sznitman (see section [T]) because we are considering several popu- 
lations with the effect that the analysis is significantly more complicated. Our 
hypotheses are also more general, e.g., the drift and diffusion functions are non- 
trivial and satisfy the general condition |(H4)| which is more general than the 
usual linear growth condition. Also they are only assumed locally (and not glob- 
ally ) Lipschitz continuous to be able to deal, for example, with the Fitzhugh- 
Nagumo model. A locally Lipschitz continuous case was recently addressed in a 
different context for a model of swarming in [7] . Proofs of our results for some- 
what stronger hypotheses than ours and in special cases are scattered in the 
literature. Our main contribution is that we provide a complete self-sufficient 
proof in a fairly general case by gathering all the ingredients that are required for 
our neuroscience applications. In particular, the case of the Fitzhugh-Nagumo 
model where the drift function does not satisfy the linear growth condition in- 
volves a generalization of previous works by using a Has'minskii-like argument 
based on the presence of a Lyapunov function to prove existence of solutions 
and non-explosion. 

The simulation of these equations can itself be very costly. We hence ad- 
dressed in section [4] numerical methods to compute the solutions of these equa- 
tions, in the probabilistic framework using the convergence result and standard 
integration methods of differential equations, or in the Fokker-Planck frame- 
work. The simulations performed for different values of the external input cur- 
rent parameter and one of the parameters controlling the noise allowed us to 
show that the spatio-temporal shape of the probability density function describ- 
ing the solution of the McKean- Vlasov- Fokker-Planck equation was sensitive to 
the variations of these parameters, as shown in figures [7| and [Sj However, we 
did not address the full characterization of the dynamics of the solutions in the 
present article. This appears to be a complex question that will be the subject 
of future work. It is known that for different McKean- Vlasov equations sta- 
tionary solutions of these equations do not necessarily exist and, when they do, 
are not necessarily unique (see [40]). A very particular case of these equations 
was treated in [76] where the authors consider that the function fa is linear, 
constant and ba(3{x,y) = Si3{y). This model, known as the firing-rate model, is 
shown in that paper to have Gaussian solutions when the initial data is Gaus- 
sian, and the dynamics of the solutions can be exactly reduced to a set of 2P 
coupled ordinary differential equations governing the mean and the standard 
deviation of the solution. Under these assumptions a complete study of the 
solutions is possible, and the dependence upon the parameters can be under- 
stood through bifurcation analysis. The authors show that intrinsic noise levels 
govern the dynamics, creating or destroying fixed points and periodic orbits. 

The mean-field description has also deep theoretical implications in neuro- 
science. Indeed, it points towards the fact that neurons encode their responses 
to stimuli through probability distributions. This type of coding was evoked by 
several authors [58 , and the mean- field approach shows that under some mild 
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conditions this phenomenon arises: ah neurons belonging to a particular popu- 
lation can be seen as independent realizations of the same process, governed by 
the mean-field equation. The relevance of this phenomenon is reinforced by the 
fact that it has recently been observed experimentally that neurons had corre- 
lation levels significantly below what had been previously reported [26]. This 
independence has deep implications on the efficiency of neural coding which 
the propagation of chaos theory accounts for. To illustrate this phenomenon 
we have performed the following simulations. Considering a network of 2, 10 
and 100 Fitzhugh-Nagumo neurons, we have simulated 2000 times the network 
equations over some time interval [0, 100]. We have picked at random a pair of 
neurons and computed the time variation of the cross-correlation of the values 
of their state variables. The results are shown in figure [9j It appears that the 
propagation of chaos is observable for relatively small values of the number of 
neurons in the network, thus indicating once more that the theory developed in 
this paper in the limit case of an infinite number of neurons is quite robust to 
finite size effect^ 




Time 



Figure 9: Variations over time of the cross-correlation of the (V^w^y) variables 
of two Fitzhugh-Nagumo neurons in a network. Top left: 2 neurons. Top right: 
10 neurons. Bottom: 100 neurons. The cross-correlation decreases steadily with 
the number of neurons in the network. 

The present study develops theoretical arguments to derive the mean-field 
equations resulting from the activity of large neurons ensembles. However, the 
rigorous and formal approach developed here does not allow direct characteri- 

^Note that we did not estimate the correlation within larger networks since, as predicted 
by theorem^ it will be smaller and smaller requiring an increasingly large number of Monte 
Carlo simulations. 
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zation of brain states. The paper however opens the way to rigorous analysis of 
the dynamics of large neuron ensembles through derivations of different quanti- 
ties that may be relevant. A first approach could be to derive the equations of 
the successive moments of the solutions. Truncating this expansion would yield 
systems of ordinary differential equations that can give approximate informa- 
tion on the solution. However, the choice of the number of moments taken into 
account is still an open question that can raise several deep questions [49j. 
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A Proof of theorem |4] 

In this appendix we prove the convergence of the network equations towards 



the mean-field equations (20) and of the propagation of chaos property. The 
proof follows standard proofs in the domain as generally done in particular by 
Tanaka or Sznitman [69] [65] , adapted to our particular case where we consider 
a non-zero drift function and a time and space dependent diffusion function. It 
is based on the very powerful coupling argument, that identifies the almost sure 
limit of the process X* as the number of neuron tends to infinity, as popularized 
by Sznitman in [67], but whose idea dates back from the 70s (for instance 
Dobrushin uses it in [25]). This process is exactly the solution of the mean- 
field equation driven by the same Brownian motion as X* and with the same 
initial condition random variable. In our case, this leads us to introduce the 
sequence of independent stochastic processes (X^*)i=i...Ar having the same law 
as X", a = p{i)^ solution of the mean- field equation: 

p 

dXi = Ut, XI) dt^Yl ^z[ba^{Xi, Z2)\ dt + ga{t, Xi) dWi 

7=1 

p 

-^Y^Mz[f3a^{XlZ^)]dBl (31) 

7=1 

with initial condition Xq = Xq, the initial condition of the neuron i in the 
network, which were assumed to be independent and identically distributed. 
{W^) and (Bl) are the Brownian motions involved in the network equation 
([T8|. As previously Z = (Z^, • • • , Z^) is a process independent of X that has 
the same law. Denoting as previously by mf the probability distribution of Xf 



solution of the mean- field equation (20), the law of the collection of processes 
(Xl^) for some fixed k G N*, namely m^^*^^ (g) . . . (g) mP^'^^\ is shown to be the 



limit of the processes (X^) solution of the network equations (18) as X goes to 
infinity. 

We recall for completeness theorem [4j 
Theorem 4. Under the assumptions \ ( Hl^ to \(H4)\ the following holds true: 
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• Convergence: For each neuron i of population a, the law of the multi- 
dimensional process X*'^ converges towards the law of the solution of the 
mean-field equation related to population a, namely . 

• Propagation of chaos: For any /c G N*^ and any k-uplet {ii^ . . . ^ik), the 
law of the process (X^*^'^, . . . ,X^*'"^,t < T) converges towards m^^^^^ (g) 
... m^''^'^\ i.e. the asymptotic processes have the law of the solution of 
the mean-field equations and are all independent. 



Proof. On our way we also prove that 



max XE 

i=l---N 



sup 

s<T 



|X!'^-X!lP 



< oo, (32) 



which imphes in particular convergence in law of the process {Xl'^ ^t < T) 



towards (X^",t < T) solution of the mean- field equations (20). 

The proof basically consists in thoroughly analyzing the difference between 
the two processes as X tends to infinity. The difference is the sum of eight 
terms (we dropped the index X for the sake of simplicity of notations) denoted 
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At through Hf. 

XI - X\ = f Us, Xi) - Us, Xi) ds + f gais, Xi) - g^{s, Xl) dWl 
Jo Jo 



At 



1 ^ 

" V ' 

Ct 

P rt . 



Dt 



7=1 ^ .7 = 1 



Et 



J2 [ ^-^(^i'^i) - pc^^^uxi) dBi 

7=1 ^^, = 1 



Ft 



fwll /?"7(^:. ^i) - PM,xi) dBp 

7=1 ^^^i=l 

E / TrE/^-7W.^f)-Ez[/3«7(^:,^7)^^r 

7=1-^0 ^V7 



(33) 



It is important to note that the probability distribution of these terms does not 
depend on the neuron i. We are interested in the hmit, as N goes to infinity, of 
the quantity ]E[supg<2^ 11-^]'^ ~ ^slP]- We decompose this expression into the 
sum of the eight terms involved in equation (33) using Holder's inequality and 
upperbound each term separately. The terms At and Bt are treated exactly as 
in the proof of theorem [2] We start by assuming that / and g are uniformly 
globally K Lipschitz-continuous with respect to the second variable. The locally- 
Lipschitz case is treated in the same manner as done in the proof of theorem [2] by 
i) stopping the process at time r[/, ii) using the Lipschitz-continuity of / and g in 
the ball of radius U and iii) by a truncation argument and using the almost sure 
boundedness of the solutions extending the convergence to the locally-Lipschitz 
case. 
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As seen previously we have: 



Now for Ct 



nsnp\\Af]<K'T / n^np\\X:-X:f]ds 

s<t Jo u<s 

E[sup||i?,f] <4i^2 f\[snp\\Xi-Xif]ds 

s<t Jo u<s 



\\Cs 



E /' ]f E ^u) - ba^iXl Xi) du 

(Cauchy-Schwarz) <TP E E ll^"7(^«. ^l) " K^iK, || du 

7' 

Jo 



(assumption [(H2)t < TPL^ / ||X; - X^f dix 
Therefore 

sup IP <TPLM \\Xi-Xifds 

s<t Jo 
rt 



E 



sup \\Cs 

,s<t 



<TPL^ / E 







sup||X;-X^ 



(is. 



Similarly for Dt: 



sup||Z),||^ <T 

s<t JO 



P . _ 

7=1 ^ j=i 



ds 



rt P 



(Cauchy-Schwartz) <PT \ - 6,^(X:,Xf)f j 

rt ( P I _ \ 

(assumption p2)]) < PTL^ \Y.mY. W^s - ^if ds 

\7=1 ^.7 = 1 / 



Hence we have 
E sup||L>5 

s<t 



t / p 



^PT^'I^ lE^.E^N-^il 

° V7=l ^j=l 
rt ( P . r 

^^^^'/ E]^ E^ sup||xi-x^ 



ds 



ds 
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Therefore 



E 



sup \\Ds 

s<t 



/ m 

Jo 3=1 



max IE 

...AT 



ds. 



The terms Ft and Gt are treated in the same fashion, but instead of using the 
Cauchy- Schwartz inequahty make use of the Bur kholder-Davis- Gundy martin- 
gale moment inequahty. For F^, in detaih 



E[sup||F,|p](Cauchy-Schwartz) < 4P V / E 

s<t 7=1-^0 



^ .1 = 1 



P pt 



(Cauchy-Schwartz) < ^pfZ f ^ E ^ \\\Pa^{XlXl) - P^^{XlXl) 
(assumption p2)t < 4.P^L^ E [||. 



\X\ — Xl 



ds 



< 4L^P^ I E 
/o 



sup||x;-x 

u<s 



i ||2 



Similarly for Gt we obtain: 



E 



sup \\Gs 

s<t 



' / ma 

Jo 3=1- 



< AL^P I max E 

■N 



sup \\Xi-Xi\ 

P<u<s 



ds. 



ds 



We are left with the problem of controlling the terms Et and Ht that involve 
sums of processes with bounded second moment thanks to proposition [3] and 
assumption |(H3)| We have: 



E 











= E 


sup 


_S<t 




s<t 









.s P . 



P nt 

(Cauchy-Schwartz) < TP V / E 



3 = 1 



ds 



and using the Burkholder-Davis-Gundy martingale moment inequality 



E 



sup \\Hs 

s<t 



p pt 



<4PV / E 



- Pa,{XlXl) - Ez [Pa,{Xl Z2 

^ 3 = 1 



ds 



Each of these two expressions involves an expectation which we write 

21 



E 



- ^6(x:,xf)-Ez[e(x:,z7)] 



3 = 1 
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where G {^^7, /^a^}, and expand as 



^ E[(6(X:,Xf)-Ez [e(X:,Z7)]f (6(X:,X,^)-Ez[6(X:,Z7)])] 

All the terms of the sum corresponding to indexes j and k such that the three 
conditions j ^ k ^ i and j ^ k are satisfied are null, since in that case XL 
Xl^ X^ and are independent and with the same law for p{j) = p{k) = 
In effect, denoting by the measure of their common law, we have: 

E [(e(x:,xi) -Ez[e(x:,z2)]r(e(xi,x,^) -Ez[e(xi,z7)])] = 

n^iXlXif ( e{Xlz)m1{dz)]- 



m]{dz)e{Xi,X'^ 



nj Q{Xl,zfm1{dz) j 



m2{dz) j e{Xi,z)m]{dz)], 

expanding further and renaming the second z variable to y in the last term we 
obtain 



/// 



G{x, y f&ix, z) m2{dx)m2{dy)m2{dz)- 

1 1 &{x,yf J Q{x,z)mZ{dz)mZ{dx)mZ{dy)- 
0(x, z)^ m2{dz)Q{x^ y) m2{dx)m2{dy)-\- 

Q{x^ z)^mj{dz) / Q{x^y)m2{dy) m2{dx) 



which is indeed equal to by Fubini theorem. 

Therefore, there are no more than 3N^ non-null terms in the sum, and 
all the terms have the s ame va lue (that depends on B), which is bounded by 
lemma [3] and assumption (H3) We denote by C/3 the supremum of these 2P^ 
values for O G {ba-f^Pa-/} across all possible pairs of populations, and by N^i^ 
the smallest value of the N^^ 7 = 1 • • • P. We have shown that 



E 



Finally we have: 



sup 11^, f 


and E 


sup||i7,|p 


_S<t 




_s<t 



< 



Nrr 



max E 

i=l---N 



sup \\Xl - XI 

,s<t 



max E 



i=i---^ \-u<s 



sup||X^-X^| 



du ■ 



Nr. 



^Note that i / j and i ^ k as soon as p{i) / p{j) = p(k) = 7. In the case where p(i) = 7, 
it is easy to check that when j (respectively k) is equal to i, all terms such that k ^ j 
(respectively j ^ k) are equal to 0. 
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for some positive constants Ki and K2. Using Gronwall's inequality, we obtain: 



max IE 

i=l---N 



sup XI -XI 

s<t 



< 



Nr. 



(34) 



for some positive constant K^, term that tends to zero as TV goes to infinity 
proving the propagation of chaos property. In order to show a convergence with 
speed 1 / ^/N as stated in the theorem, we use the fact 



max 7VE 

i=l---N 



sup ii^:-^ - xi 

,s<T 



N 



and the righthand side of the inequahty is bounded for ah A^s because of the 
hypothesis hm^v^oo ^ = ^ (O5 1) for a = 1 • • • P. 
This ends the proof. 

□ 
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